rm(list=ls(all=TRUE))

path <- "" #here you need to adjust the path macro to the place on your computer
setwd(path)

library(haven)
library(dplyr)
library(ggplot2)
library(Hmisc)
library(lme4)
library(texreg)
library(effects)
library(estimatr)
library(survey)
library(sjPlot)
library(lmtest)
library(plm)
library(broom)

LevNum <- function(x){
  c(round(quantile(x, na.rm=T)[2], digits=5),
    round(mean(x, na.rm=T), digits = 5),
    round(quantile(x, na.rm=T)[4], digits=5))
}

get.or.se <- function(model) {
  broom::tidy(model) %>% 
    mutate(or = exp(estimate),
           var.diag = diag(plm::vcovHC(model, method="arellano",type="HC1",cluster="country")),
           or.se = sqrt(or^2 * var.diag),
           conf.low.or = or-1.405054*or.se, #84 percent, non overlapping ci, for 95 use 1.959947
           conf.high.or = or+1.405054*or.se) %>% 
    select(term, or, conf.low.or, conf.high.or) %>% 
    rename(estimate = or,
           conf.low = conf.low.or,
           conf.high = conf.high.or)
}

dat <- read_dta(paste0(path, "\\chapter_wide.dta"))

#################################
###### Table 1, appendix B ######
#################################

dat$iso2code <- as_factor(dat$iso2code)
dat$country <- as_factor(dat$country)

# adjusting the sample, dropping parties which have not sponsored 
# protests and cabinets under which no protest has happened.

table(dat$country)
length(unique(dat$country))


# low and high contexts

cab_cntr <-  dat %>% 
  group_by(iso2code, cabinet_name) %>%
  summarise_at(vars(aglpopw_event, wparties), funs(sum(., na.rm=TRUE)))

cab_cntr <- within(cab_cntr, prt_share <- wparties*100/aglpopw_event)
table(cab_cntr$prt_share)
cab_cntr$cabinet_name[cab_cntr$prt_share>100]
quantile(cab_cntr$prt_share)[3]

diff.cntry <- c("BE" ,"DK" ,"FI" ,"FR" ,"GR" ,"IE" ,"IT" ,"LU" ,"MT" ,"NL" ,"NO" ,"PT" ,"RO" ,"SI" ,"UK")
nodiff.cntry <- c("AT" ,"BG" ,"CY" ,"CZ" ,"EE" ,"DE" ,"HU" ,"IS" ,"LV" ,"LT" ,"PL" ,"SK" ,"ES" ,"SE" ,"CH")


cab_cntr$prt_share_cat[as.character(cab_cntr$iso2code) %in% diff.cntry] <- 0
cab_cntr$prt_share_cat[as.character(cab_cntr$iso2code) %in% nodiff.cntry] <- 1

cab_cntr$prt_share_cat <- as.factor(cab_cntr$prt_share_cat)
cab_cntr$prt_share_cat <- relevel(cab_cntr$prt_share_cat, "0")
cab_cntr$all_evts <- as.numeric(cab_cntr$aglpopw_event)

dat_upper_lev_cat <- subset(cab_cntr, select=c("iso2code",
                                               "cabinet_name",
                                               "prt_share",
                                               "prt_share_cat",
                                               "all_evts"))
dat <- merge(dat, dat_upper_lev_cat)
rm(dat_upper_lev_cat, cab_cntr)

# Event level regression

dat$economic <- as.factor(dat$economic)
dat$NSM <- as.factor(dat$NSM)
dat$ANSM <- as.factor(dat$ANSM)
dat$political <- as.factor(dat$political)
dat$unions <- as.factor(dat$unions)
dat$otherorg <- as.factor(dat$otherorg)
dat$petition[grepl("petition", dat$action)] <- 1
dat$petition[is.na(dat$petition)] <- 0
dat$petition <- as.factor(dat$petition)
dat$strike[grepl("strike", dat$action)] <- 1
dat$strike[is.na(dat$strike)] <- 0
dat$strike <- as.factor(dat$strike)
dat$confrontation[grepl("blockade", dat$action)] <- 1
dat$confrontation[is.na(dat$confrontation)] <- 0
dat$confrontation <- as.factor(dat$confrontation)
dat$violence[grepl("violent", dat$action)] <- 1
dat$violence[is.na(dat$violence)] <- 0
dat$violence <- as.factor(dat$violence)
dat$demo[grepl("demonstration", dat$action)] <- 1
dat$demo[is.na(dat$demo)] <- 0
dat$demo <- as.factor(dat$demo)
dat$aft_elect_prt <- as.factor(dat$aft_elect_prt)
dat$bef_elect_prt <- as.factor(dat$bef_elect_prt)
dat$ctr_year <- as.factor(dat$ctr_year)
dat$all_evts <- as.numeric(dat$all_evts)


m.1 <- glm(parties ~ 
             economic * prt_share_cat +
             NSM * prt_share_cat +
             ANSM * prt_share_cat +
             political * prt_share_cat +
             aglpopw_part_all* prt_share_cat +
             unions * prt_share_cat +
             otherorg * prt_share_cat +
             petition* prt_share_cat +
             strike * prt_share_cat +
             confrontation * prt_share_cat +
             violence * prt_share_cat +
             demo * prt_share_cat +
             aft_elect_prt * prt_share_cat +
             bef_elect_prt * prt_share_cat + 
             poly(all_evts, 2):prt_share_cat +
             country,
           family = "binomial",
           weights=aglpopw_event,
           data=dat)

dat_nodiff <- dat[dat$prt_share_cat==1,] #many parties creates a non-differentiated landscape
dat_diff <- dat[dat$prt_share_cat==0,]

m_diff_1 <- glm(parties ~ economic +
             NSM + ANSM + political + 
             aglpopw_part_all +
             unions + otherorg + petition +
             strike + confrontation + violence +
             demo + aft_elect_prt + bef_elect_prt +
             poly(all_evts, 2) + country,
             family = "binomial",
             weights=aglpopw_event,
             data=dat_diff)

m_nodiff_1 <- glm(parties ~ economic +
                  NSM + ANSM + political + 
                  aglpopw_part_all +
                  unions + otherorg + petition +
                  strike + confrontation + violence +
                  demo + aft_elect_prt + bef_elect_prt +
                  poly(all_evts, 2) + country,
                family = "binomial",
                weights=aglpopw_event,
                data=dat_nodiff)

m.1.cf <- coeftest(m.1, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))
m_diff_1.cf <- coeftest(m_diff_1, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country")) 
m_nodiff_1.cf <- coeftest(m_nodiff_1, vcov.=function(x) vcovHC(x, method="arellano", type="HC1", cluster="country"))

htmlreg(list(m_diff_1, m_nodiff_1, m.1),
          custom.model.names=c("High Diff.", "Low Diff.", "Full Sample"),
          override.se = list(m_diff_1.cf[,2], m_nodiff_1.cf[,2], m.1.cf[,2]), 
          override.p=list(m_diff_1.cf[,4], m_nodiff_1.cf[,4], m.1.cf[,4]), 
          omit.coef="country",
          caption = "",
          single.row=TRUE,
          custom.note="*** p < 0.001, ** p < 0.01, * p < 0.05<br/>All models include country fixed effects.<br />Standard errors clustered by countries.",
          custom.coef.map = list("(Intercept)"="Intercept",
                                 "economic1" = "Economic",
                                 "NSM1" = "Cultural lib.",
                                 "ANSM1" = "Cultural cons.",
                                 "political1" = "Political",
                                 "aglpopw_part_all" = "Participants",
                                 "unions1" = "Unions",
                                 "otherorg1" = "Other org.",
                                 "petition1" = "Petition",
                                 "strike1" = "Strike",
                                 "confrontation1" = "Confrontation",
                                 "violence1" = "Violence",
                                 "demo1" = "Demonstration",
                                 "daft_parl_cycle" = "Electoral cycle",
                                 "aft_elect_prt1" = "After election",
                                 "bef_elect_prt1" = "Before election",
                                 "poly(all_evts, 2)1"="Overall protest",
                                 "prt_share_cat0:poly(all_evts, 2)1"="Overall protest",
                                 "poly(all_evts, 2)2"="Overall protest (sq)",
                                 "prt_share_cat0:poly(all_evts, 2)2"="Overall protest (sq)",
                                 "prt_share_cat1" = "Differentiation",
                                 "economic1:prt_share_cat1" = "Economic*Diff",
                                 "prt_share_cat1:NSM1" = "Cult. lib.*Diff",
                                 "prt_share_cat1:ANSM1" = "Cult. cons.*Diff",
                                 "prt_share_cat1:political1" = "Political*Diff",
                                 "prt_share_cat1:aglpopw_part_all" = "Participants*Diff",
                                 "prt_share_cat1:unions1" = "Unions*Diff",
                                 "prt_share_cat1:otherorg1" = "Other org.*Diff",
                                 "prt_share_cat1:petition1" = "Petition*Diff",
                                 "prt_share_cat1:strike1" = "Strike*Diff",
                                 "prt_share_cat1:confrontation1" = "Confrontation*Diff",
                                 "prt_share_cat1:violence1" = "Violence*Diff",
                                 "prt_share_cat1:demo1" = "Demonstation*Diff",
                                 "prt_share_cat1:daft_parl_cycle" = "Elect. cycle*Diff",
                                 "prt_share_cat1:aft_elect_prt1" = "After elec.*Diff",
                                 "prt_share_cat1:bef_elect_prt1" = "Before elec.*Diff",
                                 "prt_share_cat1:poly(all_evts, 2)1" = "Overall protest * Diff.",
                                 "prt_share_cat1:poly(all_evts, 2)2" = "Overall protest (sq) * Diff."),
        file = paste0(path, "\\evt_regression.html"))
